.. _tut6_ext1:
TUTORIAL 6 : Exercise 1
=======================
:Author: Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)
Reproducing *U(r)* with EE & WL/RE methods
------------------------------------------
In the limit of *infinite dilution*, i.e. when the solute concentration (or number density) is negligible,
the potential of mean force approaches the underlying effective interaction for a pair of solute particles.
This implies that if only two particles are present in the simulation box,
.. math::
{\rm PMF}=-\ln g(r)\ \to\ \beta\sum_k U_k(r)
In this exercise we will see that it is actually the case.
Exercise 1.1: Using the EE method
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
First, we will consider a pair of LJ particles in vacuum.
Change to directory **tutorial6-1/fed_2ljp_ee/**, where you will find the input files for this exercise.
Check the **cutoff**, energy **units** and **particle diameter** in FIELD file:
.. code-block:: html
:linenos:
Lennard-Jones - 2 particles
CUTOFF 12.0
UNITS K
NCONFIGS 1
atoms 1
Lj core 1.0 0.00
MOLTYPES 1
test
MAXATOM 2
vdw 1
Lj core Lj core lj 1.000000 3.000
FINISH
CLOSE
Check the **FED input** in CONTROL file:
.. code-block:: html
:linenos:
# FED block is started by 'use fed ' (must be in 'use' section on top)
use fed generic 1
# FED [ ]
fed method EE 0.5 1.125 500000 #3 10 250 #0.333
# FED order [param]
fed order parameter com2 250 2.5 27.5 1
com1 mol 1 atom 1 # use COM1 for a group of molecules or atoms within molecules
com2 mol 1 atom 2 # use COM2 for a group of molecules or atoms within molecules
com sampling correction 1 # entropy on a spherical surface (on the fly)
#com sampling correction 2 # entropy in a spherical bin volume (a posteriori)
fed order parameter done
# close FED block started by 'use fed' directive
fed done
**NOTE:** Since at every MC step one of the FED-involved particles is moved,
the stride for FED attempts can only be '1' (it will default to that otherwise)!
We start an EE iteration with 500,000 MC steps per EE update cycle initially. After two first cycles
the update stride will be doubled every time the FED bias is updated (and stored in FEDDAT._
file), i.e. we are doing a sequence of 2 x 500,000; then 1,000,000; then 2,000,000; and so on...
The EE method requires a **damping factor**, :math:`\eta`, for the energy feedback term, :math:`\Delta U_{bias}(R)`
(see above), which in this case is :math:`\eta_0=0.5`. The damping factor is scaled up after every EE iteration
by the **scaling coefficient** as :math:`\eta_{k+1}=1.125\ \eta_k` which is done with the same frequency as defined
by the main update stride.
The subsection defining the order parameter is started by directive '**fed order** ...' and closed
by '**fed order [param] done**'. The first (opening) record in this subsection specifies: the parameter
grid number, the allowed range of variation (working window), and the power of incremental bin size
(if greater than 1 in absolute value).
The specs for the **two COM groups** are self-evident: each contains a single atom (indeces 1 and 2)
in the same 'molecule 1' (under name 'test' in FIELD and CONFIG). The third directive in the 'fed order'
subsection specifies how the bias unfolding correction is to be accounted for. In this case we opted
for the 'on the fly' correction that accounts for the *ideal entropy on a sphere* of the actual (current)
radius, :math:`\sim 2\ln R`.
Further in CONTROL file the number of MC steps and strides for collecting stats should be commensurate
with those for FED:
.. code-block:: html
:linenos:
temperature 1.0 # Reduced temperature
nocoul all # No electrostatics
#distewald # Distributed Ewald sums in parallel runs
equilibration 0 # Equilibration period (not for FED calcs)
steps 4000000 # Total number of MC steps 2^(-1)*
print 1000000 # Frequency of printing to OUTPUT
stats 100000 # frequency of printing stats to PFILE
check 1000000 # Check rolling energies vs recalculated from scratch
sample rdfs 300 30.0 100 # Sampling grid, bin size & frequency for RDF calculation
sample coords 10000 # Frequency of storing configuration in trajectory
acceptatmmoveupdate 1000000000 # DO NOT update max atom displacement
maxatmdist 10.000 # Max atom displacement
move atoms 1 100 # Move atoms
Lj core # Name & type of atom to move
start simulation
Run the simulation (serial version) and check the output::
[tutorial6-1/fed_2ljp_ee]$ ls FED*
FEDDAT.000_001 FEDDAT.000_002 FEDDAT.000_003 FEDDAT.000_004
Launch gnuplot and plot the FED data, which is the second column in FEDDAT files (to quit gnuplot type 'q')::
[tutorial6-1/fed_2ljp_ee]$ gnuplot
gnuplot> plot 'FEDDAT.000_001' u 1:2 w l t "EE iter 1", 'FEDDAT.000_002' u 1:2 t "EE iter 2" w l, \
'FEDDAT.000_003' u 1:2 w l t "EE iter 1", 'FEDDAT.000_004' u 1:2 w l t "EE iter 4"
Now let's add the known target (LJ potential) and zoom in by providing X and Y ranges::
gnuplot> lj(arg) = 4*( (3.0/arg)**12 - (3.0/arg)**6 )
gnuplot> plot [x=2:10] [y=-1.5:1.5] 'FEDDAT.000_001' u 1:2 w l t "EE iter 1", \
'FEDDAT.000_002' u 1:2 t "EE iter 2" w l, 'FEDDAT.000_003' u 1:2 w l t "EE iter 3", \
'FEDDAT.000_004' u 1:2 w l t "EE iter 4", 'FEDDAT.000_001' u 1:(lj($1)) w l lc 'black' lt 0 t "Target U(r)"
.. figure:: Images-FED/FED-2ljp-ee-iters4.png
:width: 640px
Finally, also check the EE population histograms, i.e. probability distributions, which is the third column in FEDDAT files::
gnuplot> plot [] [y=0:12000] 'FEDDAT.000_001' u 1:3 w l t "EE iter 1", \
'FEDDAT.000_002' u 1:3 t "EE iter 2" w l, 'FEDDAT.000_003' u 1:3 w l t "EE iter 3", \
'FEDDAT.000_004' u 1:3 w l t "EE iter 4"
.. figure:: Images-FED/FED-2ljp-ee-iters4-probs.png
:width: 640px
What might go wrong?
^^^^^^^^^^^^^^^^^^^^
Questions to answer:
* Was the entire order parameter range covered appropriately?
* **Are the histograms sufficiently flat?**
* Can one improve of the stats by varying :math:`\eta_0`, the scaling-up coef, or the bias update stride?
* Would smoothing upong updating the bias help?
Refinement production run is necessary!
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Let's store away the current results::
[tutorial6-1]/fed_2ljp_ee$ mkdir equil-fed1/
[tutorial6-1]/fed_2ljp_ee$ cp CON* FIELD equil-fed1/
[tutorial6-1]/fed_2ljp_ee$ mv *.0??* equil-fed1/
Now amend CONTROL file:
.. code-block:: html
:linenos:
# FED [ ]
fed method EE 0.7 1.125 2000000 #3 10 250 #0.333
# FED order [param]
fed order parameter com2 250 2.5 27.5 -1
...
steps 8000000 # Total number of MC steps 2^(-1)*
and re-run.
Ooops...
Due to '-1' at the end of '**fed order**' directive, we are asking for a seed file **FEDDAT.000**
which we forgot to provide! Of course, we meant to use the FED bias data from the latest iteration::
[tutorial6-1]/fed_2ljp_ee$ cp equil-fed1/FEDDAT.000_004 ./FEDDAT.000
restart it once again and redo the analysis of the production run.
.. figure:: Images-FED/FED-2ljp-ee-iters7-probs.png
:width: 640px
.. figure:: Images-FED/FED-2ljp-ee-iters7.png
:width: 640px
Exercise 1.2: Using the WL method with replica-exchange in *T*
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Next, we will apply the Wang-Landau method to the same system of a pair of LJ particles.
Change to directory **tutorial6-1/fed_2ljp_wl/**, where the input files are found.
We will combine the WL FED calculation with the replica-exchange (RE) parallel tempering technique
so that FED(R) will be evaluated for four temperatures at once!
To this end, the following two lines in CONTROL file are added (in the '**use**' section below the title):
.. code-block:: html
:linenos:
use ortho # use orthorhombic PBC for speed-up (if the cell is cubic)
use repexch 4 -0.2 1000 # Replica-Exchange with 4 replicas, delta T = -0.2, sampling every 1000 MC steps
The FED input to use the WL method is as follows:
.. code-block:: html
:linenos:
# FED block is started by 'use fed ' directive
use fed generic 1
# FED [ ]
fed method WL 0.002 0.5 500000 #3 10 250 #0.333
# FED order [param]
fed order parameter com2 250 2.5 27.5 1
Make sure that reading the FED bias from input is turned off, i.e. '<**power**>' is set ot 1 in '**fed order** ...' directive.
Also, make sure that the rest of the input in CONTROL corresponds to our initial run comprising 4,000,000 MC steps in total.
Similar to EE, the WL iteration initially starts with 500,000 MC steps per WL update cycle. The update stride
is doubled in the same manner as for EE: 2 x 500,000; then 1,000,000; then 2,000,000; and so on...
The WL method requires an initial **droplet** value, :math:`\delta` (recall the WL illustration), which in this case is
:math:`\delta_0=0.002`. The **scaling coefficient** acts as follows: :math:`\delta_{k+1}=0.5\ \delta_k`
with the same frequency as defined by the main update stride.
**NOTE:** Even though there are no charges (no electrostatic interactions) in this case, one has to invoke
the '**distewald**' directive for *any parallel run*!
.. code-block:: html
:linenos:
temperature 1.0 # Reduced temperature
nocoul all # No electrostatics
distewald # Distributed Ewald sums in parallel runs
equilibration 0 # Equilibration period (not for FED calcs)
steps 8000000 # Total number of MC steps 2^(-1)*
print 1000000 # Frequency of printing to OUTPUT
stats 100000 # frequency of printing stats to PFILE
check 1000000 # Check rolling energies vs recalculated from scratch
sample rdfs 300 30.0 100 # Sampling grid, bin size & frequency for RDF calculation
sample coords 10000 # Frequency of storing configuration in trajectory
acceptatmmoveupdate 1000000000 # DO NOT update max atom displacement
maxatmdist 10.000 # Max atom displacement
move atoms 1 100 # Move atoms
Lj core # Name & type of atom to move
start simulation
.. If you copied and edited the correct CONTROL file from **tutorial6-1/fed_2ljp_ee/equil-fed1/**, everything should be ready by now.
To run the replica exchange parallel job, you have to ensure that you are not running an interactive session, but instead submit the job to a parallel queue::
[tutorial6-1/fed_2ljp_wl]$ sbatch parallel.sub
Run the simulation and check the output (order can be different)::
[tutorial6-1/fed_2ljp_wl]$ ls FED*
FEDDAT.000_001 FEDDAT.000_002 FEDDAT.000_003 FEDDAT.000_004 FEDDAT.000_005
FEDDAT.001_001 FEDDAT.001_002 FEDDAT.001_003 FEDDAT.001_004 FEDDAT.001_005
FEDDAT.002_001 FEDDAT.002_002 FEDDAT.002_003 FEDDAT.002_004 FEDDAT.002_005
FEDDAT.003_001 FEDDAT.003_002 FEDDAT.003_003 FEDDAT.003_004 FEDDAT.003_005
Launch gnuplot and plot the FED data::
[tutorial6-1/fed_2ljp_wl]$ gnuplot
gnuplot> plot [x=2:10] [y=-3:2] 'FEDDAT.000_005' u 1:2 w l t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:2 t "Replica 2, WL iter 5" w l, 'FEDDAT.002_005' u 1:2 w l t "Replica 3, WL iter 5", \
'FEDDAT.003_005' u 1:2 w l t "Replica 4, WL iter 5"
.. figure:: Images-FED/FED-2ljp-wl-iters5.png
:width: 640px
Check the histograms too::
[tutorial6-1/fed_2ljp_wl]$ gnuplot
gnuplot> plot [] [y=14000:18000] 'FEDDAT.000_005' u 1:3 w l t "Replica 1, WL iter 5", \
'FEDDAT.001_005' u 1:3 t "Replica 2, WL iter 5" w l, 'FEDDAT.002_005' u 1:3 w l t "Replica 3, WL iter 5", \
'FEDDAT.003_005' u 1:3 w l t "Replica 4, WL iter 5"
.. figure:: Images-FED/FED-2ljp-wl-iters5-probs.png
:width: 640px
- **Store away the data, amend the input and run the refinement (production) simulation on your own!**
Extra exercises
^^^^^^^^^^^^^^^
- **Try to run EE/RE combination**
- **Try to add charges (+1,-1) on the LJ particles and redo the FED calculations - what do you expect to obtain?**
- **Study a more physically meaningful case:**
:ref:`tut6_ext2` - Calculating PMF for two charged colloids in ionic cloud
:ref:`tut6_ext3` - Umbrella sampling in subranges (windows) with the use of WHAM method